Pro cels_mozart_emissions_write, q, lon, lat, filename=ncName, time=time
 
  fID=ncdf_create(ncName, /clobber)
  
  LonDID  = NCDF_DIMDEF(fID, 'lon', long(n_elements(lon)))
  LatDID  = NCDF_DIMDEF(fID, 'lat', long(n_elements(lat)))
  TimeDID  = NCDF_DIMDEF(fID, 'time', /unlimited)
  
  ; Create a variable to hold the data
  lonID = NCDF_VARDEF(fId,'lon', [LonDID], /DOUBLE)
  latID = NCDF_VARDEF(fId,'lat', [LatDID], /DOUBLE)
  dateID = NCDF_VARDEF(fId,'date', [TimeDID], /LONG)
  timeID = NCDF_VARDEF(fId,'time', [TimeDID], /DOUBLE)
  qID = NCDF_VARDEF(fId, 'emissions', [LonDID, LatDID, TimeDID], /DOUBLE)
  
  ; Create some attributes (to tell about our variable)
  NCDF_ATTPUT, fId, lonId, "long_name", "longitude"
  NCDF_ATTPUT, fId, lonId, "units", "degrees_east"

  NCDF_ATTPUT, fId, latId, "long_name", "latitude"
  NCDF_ATTPUT, fId, latId, "units", "degrees_north"

  NCDF_ATTPUT, fId, qId, "long_name", "emission rate"
  NCDF_ATTPUT, fId, qId, "units", "molecules/cm2/s"
    
  NCDF_ATTPUT, fId, TimeID, "long_name", "Time"
  NCDF_ATTPUT, fId, TimeID, "units", "days since 0000-01-01 00:00:00"

  NCDF_ATTPUT, fId, DateID, "long_name", "Date"
  NCDF_ATTPUT, fId, DateID, "units", "YYYYMMDD"

  ; Leave definition mode and enter data write mode
  NCDF_CONTROL, fId, /ENDEF
  
  ; Write the data
  NCDF_VARPUT, fId, lonId, lon
  NCDF_VARPUT, fId, latId, lat
 
  NCDF_VARPUT, fId, dateID, mozdate(time)
  NCDF_VARPUT, fId, timeID, time
  NCDF_VARPUT, fId, qId, q

  ; Done
  NCDF_CLOSE, fId

End

; docformat = 'rst'
;
;+
;
; :Purpose:
;   Calculate emissions fields for MOZART.
;   
;   Generates {input_directory}/{sim_name}/emissions/emissions_totals.inp.sav file that contains
;     the pulse magnitude required in sensitivity calculation
;   
; :Inputs:
;   sim_name: (string) Simulation name
;
; :Requires:
;   Martin Shultz's regridding routine. Download at http://www.acd.ucar.edu/Applications/regrid_v2.tar.gz
;   
; :Outputs:
;   MOZART emissions files in {input_directory}/emissions/MOZART/
;   
;   Emissions totals and pulse magnitudes in {input_directory}/{sim_name}/emissions/emissions_totals.inp.sav
;   
; :History:
;   Written by: Matt Rigby, MIT, Aug 19, 2011
;
;-
pro cels_mozart_emissions, sim_name

  compile_opt idl2
  on_error, 2
  
  print, "CELS_MOZART_EMISSIONS: calculating EUM emissions and pulse magnitudes"
  
  avogadro=6.0221414e23

  pollutant=cels_get_parameter(sim_name, 'POLLUTANT')
  pollutant_short=cels_get_parameter(sim_name, 'POLLUTANT_SHORT')
  stations=cels_get_parameter(sim_name, 'STATIONS')
  StartY_IC=cels_get_parameter(sim_name, 'STARTY_IC')
  StartY=cels_get_parameter(sim_name, 'STARTY')
  EndY=cels_get_parameter(sim_name, 'ENDY')
  molmass=cels_get_parameter(sim_name, 'MOLMASS')
  ic_name=cels_get_parameter(sim_name, 'IC_NAME')
  NLR_name=cels_get_parameter(sim_name, 'NLR_NAME')

  if StartY_IC eq !null then StartY_IC=StartY

  restore, cels_filestr(/Input, sim_name + '/regions.sav')
  restore, cels_filestr(/Input, sim_name + '/state.sav')

  pulse_magnitude_ic=fltarr(N_IC)
  pulse_factor_ic=1. ;MLR
  Pulse_magnitude_ic_ID=IC_NAME

  pulse_magnitude_NLR=fltarr(N_NLRT)
  pulse_factor_NLR=0.1  ;MLR
  Pulse_magnitude_NLR_ID=StrArr(N_NLRT)

  pulse_magnitude_LR=fltarr(n_LRT)
  pulse_factor_LR=10.  ;MLR
  Pulse_magnitude_LR_ID=StrArr(N_LRT)

  ;Reference file
  ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  
  restore, cels_filestr(/Input, sim_name + '/EUM/grid.inp.sav') ;kg/m2/s
  restore, cels_filestr(/Input, sim_name + '/emissions/emissions.inp.sav') ;kg/m2/s

  ;Convert to molecules/cm2/s
  q=q/100./100.*1000./molmass*avogadro

  wh=where(q_time ge julday(1, 1, StartY_IC, 0) and q_time lt julday(1, 1, EndY, 0), count)
  if count ne (EndY - StartY_IC) then begin  ;Annual only
    message, "emissions.inp.sav' does not span entire date range"
  Endif
  q_time=q_time[wh]
  q=q[*, *, wh]

  ;Martin Shultz's regridding routine. 
  ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

  ;Define old grid
  nlon = n_elements(q_lon)
  incr = (max(q_lon) - min(q_lon))/Double(nlon-1)
  ms_lon = Obj_New('MGS_LonVector', $
                start=min(q_lon),incr=incr,nval=nlon, $
                type='regular')
  nlat = n_elements(q_lat)
  incr = (max(q_lat) - min(q_lat))/Double(nlat-1)
  ms_lat = Obj_New('MGS_LatVector', $
                start=min(q_lat),incr=incr,nval=nlat, $
                type='regular')
  o = Obj_New('MGS_RGrid', name='Source Grid', $
              Longitude=ms_lon, Latitude=ms_lat)
  if total(o->getlon() - q_lon) gt 0.1 then stop
  if total(o->getlat() - q_lat) gt 0.1 then stop
  
  ;Define new grid
  LonSize = n_elements(lon)
  lon_incr = (max(lon) - min(lon))/Double(LonSize-1)
  ms_lon = Obj_New('MGS_LonVector', $
                start=min(lon),incr=lon_incr,nval=LonSize, $
                type='regular')
  LatSize = n_elements(lat)
  lat_incr = (max(lat) - min(lat))/Double(LatSize-1)
  ms_lat = Obj_New('MGS_LatVector', $
                start=min(lat),incr=lat_incr,nval=LatSize, $
                type='regular')
  n = Obj_New('MGS_RGrid', name='Target Grid', $
              Longitude=ms_lon, Latitude=ms_lat)
  
  if total(n->getlon() - lon) gt 0.1 then stop
  if total(n->getlat() - lat) gt 0.1 then stop
  
  ;Regrid
  q=regrid(q, oldgrid=o, newgrid=n, /use_grid, /area_weighted) ;kg/m2/s
  
  ;Clean
  Obj_Destroy, o
  Obj_Destroy, n

  area=areagrid(lon, lat)

  ;Annual only:
  ; Create emissions timesteps on the first and last day of each year
  ; so that MOZART interpolates constant emissions throughout each year.
  ; Also include a time step just before the first simulation time and after the last time

  TimeSize=(EndY - StartY_IC)*2 + 2  
  Time_out=dblarr(TimeSize)
  q_out=fltarr(LonSize, LatSize, TimeSize)
  
  Time_out[0]=julday(1, 1, StartY_IC-1, 0)  ;One year before
  Q_out[*, *, 0]=q[*, *, 0]
  for yi=StartY_IC, EndY-1 do begin
    Time_out[(yi - StartY_IC)*2+1]=julday(1, 1, yi, 0)
    Time_out[(yi - StartY_IC)*2+2]=julday(12, 31, yi, 0)
    Q_out[*, *, (yi - StartY_IC)*2+1]=q[*, *, yi - StartY_IC]
    Q_out[*, *, (yi - StartY_IC)*2+2]=q[*, *, yi - StartY_IC]
  Endfor
  Time_out[-1]=julday(12, 31, EndY, 0)    ;One year later
  Q_out[*, *, -1]=q[*, *, EndY - StartY_IC - 1]
  Q=temporary(q_out)
  time=temporary(time_out)
  date=mozdate(time)
  TimeSize=n_elements(time)
  
  caldat, time, dummy, dummy1, year

  ;Write reference emissions files
  fname=cels_filestr(/Input, sim_name + '/emissions/MOZART/emissions.' + pollutant_short + '.surface.CELS.nc')
  cels_mozart_emissions_write, Q, lon, lat, time=time, filename=fname

  ;Write initial conditions
  ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  
  If StartY_IC lt StartY then begin
    For xi=0, n_IC-1 do begin
      Q_ic=fltarr(LonSize, LatSize, TimeSize)
      Case strmid(ic_name[xi], 0, 3) of
        'IC': begin
          pulse_magnitude_ic[xi]=1. ;This will be used in MOZART IC routine
          wh=0
        end
        'ICQ': begin
          wh=where(year eq StartY_IC + (fix(strmid(ic_name[xi], 3, 1))-1), count)
          if count ne 0 then begin
            Q_IC[*, *, wh]=Q[*, *, wh]*pulse_factor_ic
            pulse_magnitude_ic[xi]=$
              total(q_ic[*, *, wh[0]]*area)*100.*100./avogadro*molmass/1000.  ;kg/s
          endif else begin
            message, StrCompress('Missing ' + ic_name[xi])
          endelse
        end
        else: pulse_magnitude_ic[xi]=1.
      endcase
      fname=cels_filestr(/Input, sim_name + '/emissions/MOZART/emissions.' + $
        pollutant_short + ic_name[xi] + '.surface.CELS.nc')
      cels_mozart_emissions_write, Q_IC, lon, lat, time=time, filename=fname
    endfor
  endif else begin
    for xi=0, n_IC-1 do begin
      pulse_magnitude_ic[*]=1. ;This will be used in MOZART IC routine      
      fname=cels_filestr(/Input, sim_name + '/emissions/MOZART/emissions.' + $
        pollutant_short + ic_name[xi] + '.surface.CELS.nc')
      cels_mozart_emissions_write, Q, lon, lat, time=time, filename=fname
    endfor
  endelse
   
  ;Local regions
  ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  
  lon_domain=lon_domain_EUM
  lat_domain=lat_domain_EUM
  
  ;Create map of local regions
  regions=intarr(n_elements(lon), n_elements(lat))
  for lri=0, n_LR-1 do begin
    if lon_domain[0, lri] gt lon_domain[1, lri] then begin
      whlon=where(lon gt lon_domain[0, lri] or lon lt lon_domain[1, lri])
      whlat=where(lat gt lat_domain[0, lri] and lat lt lat_domain[1, lri])
      for lonI=0, n_elements(lon)-1 do for latI=0, n_elements(lat)-1 do begin
        if total(whlon eq lonI)*total(whlat eq latI) gt 0. then regions[lonI, LatI]=lri+1
      endfor
    endif else begin
      whlon=where(lon gt lon_domain[0, lri] and lon lt lon_domain[1, lri])
      whlat=where(lat gt lat_domain[0, lri] and lat lt lat_domain[1, lri])
      for lonI=0, n_elements(lon)-1 do for latI=0, n_elements(lat)-1 do begin
        if total(whlon eq lonI)*total(whlat eq latI) gt 0. then regions[lonI, LatI]=lri+1
      endfor
    endelse
  endfor

  ;Annual perturbed maps
  For yi=StartY, EndY-1 do begin  ;Annual only
    for lri=0, n_LR-1 do begin
      q_region=fltarr(LonSize, LatSize, TimeSize)

      wh=where(year eq yi)
      for whi=0, n_elements(wh)-1 do begin
        q_ti=reform(q[*, *, wh[whi]])
        q_temp=fltarr(LonSize, LatSize)
        q_temp[where(regions eq lri+1)]=q_ti[where(regions eq lri+1)]*pulse_factor_lr
        q_region[*, *, wh[whi]]=q_temp
      endfor
  
      if yi ge 2000 then year_str=string(yi - 2000, format='(I02)') else year_str=string(yi - 1900, format='(I02)')
      pulse_magnitude_lr_ID[n_LR*(yi - StartY) + lri]=strcompress(string(lri, format='(I02)') + year_str, /remove_all)
      
      
      cels_mozart_emissions_write, q_region, lon, lat, time=time, $
        filename=cels_filestr(/Input, sim_name + '/emissions/MOZART/' + $
        'emissions.' + pollutant_short + pulse_magnitude_lr_ID[n_LR*(yi - StartY) + lri] + '.surface.CELS.nc')
      cels_mozart_emissions_write, q_region, lon, lat, time=time, $
        filename=cels_filestr(/Input, sim_name + '/emissions/MOZART/' + $
        'emissions.' + pollutant_short + 'P' + pulse_magnitude_lr_ID[n_LR*(yi - StartY) + lri] + '.surface.CELS.nc')

      pulse_magnitude_lr[n_LR*(yi - StartY) + lri]=total(q_region[*, *, wh[0]]*area)*100.*100./avogadro*molmass/1000. ;kg/s
    endfor
  endfor


  ;Non-local emissions
  ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  
  region_NLR=intarr(lonsize, latsize)
  case N_NLR of
    2: begin
      ;HEMISPHERIC
      region_NLR[*, where(lat ge 0.)]=1
      region_NLR[*, where(lat lt 0.)]=2
    end
    6: begin
      ;CONTINENTAL  ['NA', 'SA', 'EU', 'AF', 'AS', 'OC']  ;MAKE SURE THESE ARE CONSISTENT:
      CONTINENT_NAME=['NORTH AMERICA', 'SOUTH AMERICA', 'EUROPE', 'AFRICA', 'ASIA', 'OCEANIA']
      continent_i=find_country(lon, lat, name=name, /ocean, /continent)
      for ri=0, 5 do begin
        if continent_name[ri] eq 'SOUTH AMERICA' then begin
          whC=where(name eq 'CENTRAL AMERICA')
          whS=where(name eq 'SOUTH AMERICA')
          region_NLR[where(continent_i eq whC[0])]=ri+1   
          region_NLR[where(continent_i eq whS[0])]=ri+1
        endif else begin
          wh=where(name eq continent_name[ri])
          region_NLR[where(continent_i eq wh[0])]=ri+1
        endelse
      endfor
    end  
  endcase

  region_NLR[where(regions ge 1)]=0

  For yi=StartY, EndY-1 do begin  ;Annual only
    for ri=1, N_NLR do begin
      q_NLR=fltarr(LonSize, LatSize, TimeSize)

      wh=where(year eq yi)
      for whi=0, n_elements(wh)-1 do begin
        q_temp=fltarr(LonSize, LatSize)
        q_ti=reform(q[*, *, wh[whi]])
        q_temp[where(region_NLR eq ri)]=q_ti[where(region_NLR eq ri)]*pulse_factor_nlr
        q_NLR[*, *, wh[whi]]=q_temp
      endfor

      if yi ge 2000 then year_str=string(yi - 2000, format='(I02)') else year_str=string(yi - 1900, format='(I02)')
      pulse_magnitude_NLR_ID[n_NLR*(yi - StartY) + ri-1]=strcompress(NLR_name[ri-1] + year_str, /remove_all)

      cels_mozart_emissions_write, q_NLR, lon, lat, time=time, $
        filename=cels_filestr(/input, sim_name +  $
          '/emissions/MOZART/emissions.' + pollutant_short + pulse_magnitude_NLR_ID[n_NLR*(yi - StartY) + ri-1] +  '.surface.CELS.nc')
      
      pulse_magnitude_NLR[n_NLR*(yi - StartY) + ri-1]=total(q_NLR[*, *, wh[0]]*area)*100.*100./avogadro*molmass/1000. ;kg/s
    endfor
  Endfor

  ;Calculate emissions
  ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  emissions_ic=pulse_magnitude_ic/pulse_factor_ic
  emissions_NLR=pulse_magnitude_NLR/pulse_factor_NLR
  emissions_LR=pulse_magnitude_LR/pulse_factor_lr

  print, 'CELS_MOZART_EMISSIONS: These annual emissions should agree:'
  for yi=StartY, EndY-1 do begin
    wh=where(year eq yi)
    print, yi, total(q[*, *, wh[0]]*area)*100.*100./avogadro*molmass/1000.*3600.*24.*365.25/1.e6, $
      total([emissions_NLR[N_NLR*(yi - StartY):N_NLR*(yi - StartY+1)-1], $
        emissions_LR[N_LR*(yi - StartY):N_LR*(yi - StartY+1)-1]])*3600.*24.*365.25/1.e6, ' Gg/yr'
  endfor
  
  save, filename=cels_filestr(/Input, sim_name + '/emissions/emissions_totals.inp.sav'), $
    pulse_magnitude_ic, pulse_magnitude_NLR, pulse_magnitude_LR, $
    emissions_ic, emissions_NLR, emissions_LR, $
    pulse_magnitude_IC_ID, pulse_magnitude_NLR_ID, pulse_magnitude_LR_ID

end
